
do "E:/ReplicateBuild/02_code/00_environment/00_set_environment.do"

local construct_temp_build = 1
	local unshrunken_va = 1
local structural_param = 1
local va_program = 1

********Gen teacher by year VA measures 
****Merge with master build file and test accomodations file
if `construct_temp_build' == 1{
	*Get reading student-teacher match
	use "$basedata/course_mem_clean_rd.dta", clear

	keep birthdt schlcode year lea statecourse section meetingcode grade mastid teachid numstudents
	ren (*) *_rd
	 
	ren mastid_rd mastid
	ren year_rd year
	 
	tempfile temprd
	save `temprd', replace

	*Get math student-teacher match
	use "$basedata/course_mem_clean.dta", clear
	
	keep birthdt schlcode year lea statecourse section meetingcode grade mastid teachid numstudents
	ren (*) *_ma
	 
	ren mastid_ma mastid
	ren year_ma year
	 
	 
	sort mastid year
	
	*Merge reading match data
	merge 1:1 mastid year using `temprd'
	drop _m
	 
	gen birthdt = birthdt_ma
	replace birthdt = birthdt_rd if birthdt==.
	replace birthdt = . if birthdt_ma!=birthdt_rd & birthdt_ma!=. & birthdt_rd!=.
	drop birthdt_*
	 
	gen grade = grade_ma
	replace grade = grade_rd if grade==.
	replace grade = . if grade_ma!=grade_rd & grade_ma!=. & grade_rd!=.
	drop grade_*
	 
	 
	merge 1:1 mastid year using "$basedata/mb_clean.dta"
	ren _merge mrg_mb

	* keep only the standard tests
	replace score_ma = . if !inlist(ma_test_id,"RG","MA03","MA04","MA05","MA06","MA07","MA08")
	replace score_rd = . if !inlist(rd_test_id,"RG","RD03","RD04","RD05","RD06","RD07","RD08")

	merge 1:1 mastid year using "$basedata/test_accom.dta"
	ren _merge mrg_accom
	xtset mastid year
	
		* merge in behavioral outcomes
	cap drop _merge
	sort mastid year
	merge 1:1 mastid year using "$basedata/suspensions"
	drop if _m==2
	bys year: egen max_m = max(_m)
	replace anyISS = 0 if _m==1 & max_m==3
	replace anyOSS = 0 if _m==1 & max_m==3
	drop _m max_m
	
	sort mastid year
	merge 1:1 mastid year using "$basedata/attendance_mb"
	drop if _m==2
	drop _m
	
	
	******Generate age variable
	gen age =  (mdy(6,30,year)-birthdt)/365.25
	replace age=0 if age==.
	 	xtset mastid year

	egen teach_year_ma=group(teachid_ma year)
	egen teach_year_rd=group(teachid_rd year)
	egen teach_schl_ma=group(teachid_ma schlcode_ma lea_ma)
	egen teach_schl_rd=group(teachid_rd schlcode_rd lea_rd)
	gen l_grade=l.grade
	* GEN lagged scores, lagged squares, and lagged cubes
	* CREATE INTERACTIONS OF PRIOR TEST SCORES WITH GRADE LEVEL
	foreach i in ma rd {
		cap gen l_`i'=l.score_`i'
		cap gen l_`i'_sq= l_`i'^2
		cap gen l_`i'_cb= l_`i'^3
		cap gen l_`i'_miss=(l_`i'==.)
		cap xi i.l_grade|l_`i' i.l_grade|l_`i'_sq i.l_grade|l_`i'_cb  , prefix(IL`i')
	}
	***** Generate pca
		xtset mastid year

	gen log_absence_rate = log(1-attendance_rate)
	pca anyISS anyOSS log_absence_rate
	predict pc1, score
	replace pc1 = -pc1 // higher values better
	
	* GEN lagged scores, lagged squares, and lagged cubes
	* CREATE INTERACTIONS OF PRIOR TEST SCORES WITH GRADE LEVEL

	cap gen l_pc1 = l.pc1
	cap gen l_pc1_sq = l_pc1^2
	cap gen l_pc1_cb = l_pc1^3
	cap gen l_pc1_miss = (l_pc1==.)
	cap xi i.l_grade|l_pc1 i.l_grade|l_pc1_sq i.l_grade|l_pc1_cb, prefix(ILp)
	
	
	
	**Focus on grade 4-8
	drop if grade<4
	drop if grade>8
	****Generate teacher-by-school-by-studenttype-VA, teacher-by-year-VA, teacher-by-school VA
	gen lag_ma=l_ma
	gen lag_rd=l_rd


	egen teach_DISADV_yr_ma=group(teachid_ma DISADV year)
	egen teach_DISADV_yr_rd=group(teachid_rd DISADV year)
	egen teach_DISADV_schl_ma=group(teachid_ma DISADV schlcode_ma lea_ma)
	egen teach_DISADV_schl_rd=group(teachid_rd DISADV schlcode_rd lea_rd)

	egen section_size_ma=sum(1), by(year schlcode_ma lea_ma  section_ma teachid_ma)
	egen section_size_rd=sum(1), by(year schlcode_rd lea_rd  section_rd teachid_rd)
	****class demographics
		foreach ss in ma rd {
	gen achHi_`ss'=(l_`ss'>0)
	}
	gen white=(ethnic==5)
	gen black=(ethnic==4)
	gen hispanic=(ethnic==3)

	foreach ss in ma rd {
		foreach i in white black hispanic female disability DISADV ENG_LEARN {
			egen class_mn_`i'_`ss'=mean(`i'), by(year schlcode_`ss' lea_`ss'  section_`ss' teachid_`ss')
				egen schl_mn_`i'_`ss'=mean(`i'), by(year schlcode_`ss' lea_`ss')
		}
	}

	****Generate indicators for data missingness
	foreach i in grade ethnic female gifted disability MIGRANT ENG_LEARN DISADV DISADV_persist year read_accom math_accom numstudents_ma numstudents_rd {
		cap gen `i'_miss=(`i'==.)
		replace `i'=0 if `i'==.
	}
	replace numstudents_ma=section_size_ma if numstudents_ma==0
	replace numstudents_rd=section_size_rd if numstudents_rd==0
	 
	egen sidx_ma = group(schlcode_ma lea_ma)
	egen sidx_rd = group(schlcode_rd lea_rd)	 
	  
	local student_cov2 i.grade i.ethnic i.female i.gifted i.disability i.MIGRANT i.ENG_LEARN  i.DISADV_persist read_accom math_accom  *_miss i.year#i.grade age

	cap egen schid_ma= group(schlcode_ma lea_ma)
	cap egen schid_rd= group(schlcode_rd lea_rd)

	tempfile tempbuild
	save `tempbuild', replace

	if `unshrunken_va' == 1 {
	foreach ss in ma rd {
		use `tempbuild', clear
		
		local student_cov2 i.grade i.ethnic i.female i.gifted i.disability i.MIGRANT i.ENG_LEARN  i.DISADV_persist read_accom math_accom  *_miss i.year#i.grade age
		
		***Focus sample 
		qui areg pc1 IL* `student_cov2' numstudents_`ss' , absorb(teach_year_`ss')
		predict VA_1yr, d
		cap gen samp=e(sample)
		drop if samp!=1
		cap egen stu_count =sum(samp), by (teachid_`ss' year)
		capture egen tag_tch_sch_yr=tag(teachid_`ss' year schlcode_`ss' lea_`ss')
		cap egen tchr_schls=sum(tag_tch_sch_yr), by(teachid_`ss' year)
		****Drop teachers who teach SPED and ELL

		drop if class_mn_disability_`ss'>.5
		drop if class_mn_ENG_LEARN_`ss'>.5


		 
		********************************************************************
		****Generate teacher school EB VA measures

		capture egen tid=group(teachid_`ss' DISADV)
		cap egen tidy=group(year schlcode_`ss' lea_`ss' section_`ss' DISADV teachid_`ss' )
		cap egen sch_by_yr=group(schlcode_`ss' lea_`ss' year)
		cap gen ncerdc_id=teachid_`ss'
		cap gen sy=year
		merge m:1 ncerdc_id sy using "$basedata/ncerdc_experience.dta" 
		drop if _merge==2
		drop _merge


		local student_cov2 i.grade i.ethnic i.female i.gifted i.disability i.MIGRANT i.ENG_LEARN  i.DISADV_persist read_accom math_accom  *_miss i.year#i.grade age

		gen vijt=.
		gen sample=.
		foreach i in 0 1 {
			areg pc1 IL* `student_cov2' numstudents_`ss' if DISADV==`i', absorb(tidy)
			replace sample=1 if e(sample)==1
			predict temp if DISADV==`i' & e(sample)==1, d 
			predict temp1 if DISADV==`i' & e(sample)==1, r
			replace vijt=temp+temp1  if DISADV==`i' & e(sample)==1
			drop temp1 temp
		}

		save "$temp/student_level_`ss'_w_unshrunken_VA_behavioral", replace
		
	}
}

	
}

if `structural_param' == 1{
	foreach ss in ma rd { 

		use "$temp/student_level_`ss'_w_unshrunken_VA_behavioral", clear

		cap drop lea schlcode

		gen expdum = tchr_exp_pay_level
		replace expdum = 6 if tchr_exp_pay_level>6 & tchr_exp_pay_level!=.
		replace expdum = 7 if tchr_exp_pay_level==.
		assert expdum<=7 & floor(expdum)==expdum

		ren schid_`ss' schid
		ren lea_`ss' lea
		ren schlcode_`ss' schlcode
		ren teachid_`ss' teachid
		ren section_`ss' section
		ren sidx_`ss' sidx

		preserve
			gen numstudentstidy = 1
			collapse (mean) exp* tchr_exp_pay_level vijt DISADV sy (sum) numstudentstidy, by(tid sch_by_yr schid tidy)


			reghdfe vijt i.expdum [fw=numstudentstidy], absorb(tch_e0_sfe=tid sfe=schid yfe=sy) res(resid_tch_sfe)


			***predict only the fitted experience profile
			predict xb,xb
			cap gen VA_tchr_exp_sfe= xb+tch_e0_sfe

			collapse (mean) xb VA_tchr_exp_sfe sfe, by(tid sch_by_yr schid tidy)
			tempfile tempVA
			save `tempVA', replace
		restore

		sort tid sch_by_yr schid tidy
		merge n:1 tid sch_by_yr schid tidy using `tempVA'
		assert _m==3
		drop _m

		* construct class variable
		egen jtidx = group(teachid year schlcode lea)
		sort jtidx section
		gen newsection = jtidx!=jtidx[_n-1] | section!=section[_n-1]
		gen secnum = 1 if jtidx!=jtidx[_n-1]
		replace secnum = secnum[_n-1]+newsection if secnum==.

		gen j = teachid // teacher index
		gen c = secnum // class index (numbered 1 to C, reset for each teacher-year)
		gen m = 1*(DISADV==0) + 2*(DISADV==1) // student types (numered 1 to M)
		gen t = year // year
		gen r = vijt - xb - sfe // student-level residual: original residual minus experience effects minus school-year FEs
		gen alphaZ = xb // estimated experience effect (from prior step), varies by teacher-year
		gen e = tchr_exp_pay_level
		gen s = sidx

		keep j c m t r alphaZ e s schlcode lea vijt sfe

		qui summ t
		global tmin = r(min)
		global tmax = r(max)

		global nt = $tmax - $tmin

		** estimate the variance of the student error, by type

		* take the student residual, minus the classroom mean, by type
		* square the residual and take the mean across the data (mean residual is 0, so square = variance)

		bys j c m t: egen mean_resid_ct = mean(r)
		bys j c m t: egen n_resid_ct = count(r)
		gen resid_i_sq = (r-mean_resid_ct)^2
		gen resid_i_sq_c = resid_i_sq * n_resid_ct/(n_resid_ct-1)

		*Variance of student residuals
		qui summ resid_i_sq_c if m==1
		gen sigma2_eps1 = r(mean)

		qui summ resid_i_sq_c if m==2
		gen sigma2_eps2 = r(mean)

		** collapse the student residuals to the teacher-classroom-type level 
		* calculate the mean and the number of observations

		collapse (count) n_ct=r (mean) Abar=r alphaZ sigma* e vijt sfe, by(j c m t s schlcode lea)

		** estimate the covariances and variances necessary for the structural parameters

		* cov(\bar{A}_{jmct},\bar{A}_{jmc't}) for m=1,2
		* will just keep classrooms 1 and 2. can add more for precision

		preserve

			gen Abar_c1 = Abar if c==1
			gen Abar_c2 = Abar if c==2

			collapse (mean) Abar_c1 Abar_c2, by(j m t s)

			keep if Abar_c1!=. & Abar_c2!=.

			corr Abar_c1 Abar_c2 if m==1, covariance
			local sigma2_muj1 = r(cov_12)

			corr Abar_c1 Abar_c2 if m==2, covariance
			local sigma2_muj2 = r(cov_12)

		restore

		gen sigma2_muj1 = `sigma2_muj1'
		gen sigma2_muj2 = `sigma2_muj2'

		* cov(\bar{A}_{j1ct},\bar{A}_{j2c't})

		preserve

			gen Abar_c1_m1 = Abar if c==1 & m==1
			gen Abar_c1_m2 = Abar if c==1 & m==2
			gen Abar_c2_m1 = Abar if c==2 & m==1
			gen Abar_c2_m2 = Abar if c==2 & m==2

			collapse (mean) Abar_c1_m* Abar_c2_m*, by(j t s)

			* we want to make both types of comparisons (c=1,m=1) vs. (c=2,m=2) AND (c=1,m=2) vs. (c=2,m=1)
			* expand data
			expand 2, gen(obsnum)

			gen Abar_c_m_1 = Abar_c1_m1 if obsnum==1
			replace Abar_c_m_1 = Abar_c2_m1 if obsnum==0

			gen Abar_c_m_2 = Abar_c2_m2 if obsnum==1
			replace Abar_c_m_2 = Abar_c1_m2 if obsnum==0

			corr Abar_c_m_1 Abar_c_m_2, covariance
			local sigma_muj1muj2 = r(cov_12)

		restore

		gen sigma_muj1muj2 = `sigma_muj1muj2'

		* cov(\bar{A}_{j1ct},\bar{A}_{j2ct})

		preserve

			gen Abar_m1 = Abar if m==1
			gen Abar_m2 = Abar if m==2

			collapse (mean) Abar_m1 Abar_m2, by(j c t s)

			corr Abar_m1 Abar_m2, covariance
			local cov_m1m2 = r(cov_12)

		restore

		gen sigma_theta1theta2 = `cov_m1m2' - sigma_muj1muj2

		* Var(\bar{A}_{jmct})

		forv mm=1/2 {

			qui summ Abar if m==`mm'
			gen VarAbar_m`mm' = r(Var)

			gen temp`mm' = VarAbar_m`mm' - sigma2_muj`mm' - (sigma2_eps`mm' / n) if m==`mm'

			qui summ temp`mm' 
			gen sigma2_theta`mm' = r(mean)

		}

		drop temp*


		* estimate drift parameters, by comparing across years (sometimes within, sometimes across type)
		* estimate \sigma_{\mu_1,|t-s|}, \sigma_{\mu_2,|t-s|}, \sigma_{\mu_1 \mu_2,|t-s|}

		preserve

			gen Abar_m1_wgt = Abar*n_ct if m==1 // want to calculate mean residual, weighted by number of students
			gen Abar_m2_wgt = Abar*n_ct if m==2

			gen n_ct_m1 = n_ct if m==1
			gen n_ct_m2 = n_ct if m==2

			collapse (sum) n_ct_m1 n_ct_m2 Abar_m1_wgt Abar_m2_wgt, by(j t)

			gen Abar_m1 = Abar_m1_wgt / n_ct_m1
			gen Abar_m2 = Abar_m2_wgt / n_ct_m2

			egen teacher_index = group(j)
			gen teacher_odd = mod(teacher_index,2)==1

			xtset j t

			local nt = $nt

			forv ts=1/`nt' {
				gen Abar_m1_L`ts' = l`ts'.Abar_m1
				gen Abar_m2_L`ts' = l`ts'.Abar_m2
				
				corr Abar_m1 Abar_m1_L`ts', covariance
				local cov_m1m1_`ts' = r(cov_12)
				
				corr Abar_m2 Abar_m2_L`ts', covariance
				local cov_m2m2_`ts' = r(cov_12)
				
				gen x1 = Abar_m1*(teacher_odd==1) + Abar_m2*(teacher_odd==0)
				gen x2 = Abar_m2_L`ts'*(teacher_odd==1) + Abar_m1_L`ts'*(teacher_odd==0)
				corr x1 x2, covariance
				local cov_m1m2_`ts' = r(cov_12)
				
				drop Abar_m1_L`ts' Abar_m2_L`ts' x1 x2

			}
		restore

		local nt = $nt

		forv ts=1/`nt' {
			gen sigma_m1m1_`ts' = `cov_m1m1_`ts''
			gen sigma_m2m2_`ts' = `cov_m2m2_`ts''
			gen sigma_m1m2_`ts' = `cov_m1m2_`ts''
		}

		** estimate the reduced form parameters


		* start by collapsing the data so that types are in columns, not separate rows

		forv mm=1/2 {
			gen n_ct_m`mm' = n_ct if m==`mm'
			gen Abar_m`mm' = Abar if m==`mm'

		}

		collapse (mean) n_ct_m* Abar_m* alphaZ sigma* e vijt, by(j c t s schlcode lea)

		assert s!=.
		bys j t: egen mins = min(s)
		bys j t: egen maxs = max(s)
		gen multi_school = mins!=maxs
		drop mins maxs

		sort j s multi_school
		gen teach_school = 1 if j!=j[_n-1] & multi_school!=1
		replace teach_school = teach_school[_n-1] + (s!=s[_n-1]) if teach_school==. & multi_school!=1

		sort j t s multi_school
		gen teach_school_order = 1 if j!=j[_n-1] & multi_school!=1
		replace teach_school_order = teach_school_order[_n-1] + (s!=s[_n-1]) if teach_school_order==. & multi_school!=1

		save "$temp/VAstructural_`ss'_behavioral", replace

		qui summ t
		global tmin_`ss' = r(min)
		global tmax_`ss' = r(max)

		qui summ teach_school
		global smin_`ss' = r(min)
		global smax_`ss' = r(max)

		qui summ teach_school_order
		global sminorder_`ss' = r(min)
		global smaxorder_`ss' = r(max)

		use "$temp/VAstructural_`ss'_behavioral", replace

		collapse (sum) n_ct_m1 n_ct_m2 (mean) multi_school teach_school teach_school_order, by(j s t schlcode lea)

		gen p_m1 = n_ct_m1 / (n_ct_m1 + n_ct_m2)
		gen p_m2 = n_ct_m2 / (n_ct_m1 + n_ct_m2)

		drop n_ct_m*

		save "$temp/VAschoolindices_`ss'_behavioral", replace

	}
}

if `va_program' == 1{
	*******
	cap program drop blp_va
	program define blp_va

	di "Shrinking VA for Subject `1', VA Type `2', VA Type Detail `3'"

	use "$temp/VAstructural_`1'_behavioral", clear

	gen sigma_m1m1_0 = sigma2_muj1
	gen sigma_m2m2_0 = sigma2_muj2
	gen sigma_m1m2_0 = sigma_muj1muj2

	* save experience measures to merge on at end

	preserve

	if "`2'" == "career" {
		keep if t==`3'
		collapse (mean) alphaZ e sigma*, by(j t)
	}
	if "`2'" == "loY" {
		keep if t==`3'
		collapse (mean) alphaZ e, by(j t)
	}
	if "`2'" == "preY" {
		keep if t==`3'
		collapse (mean) alphaZ e, by(j t)
	}
	if "`2'" == "loS" | "`2'" == "wiS" {
		keep if t==`3'
		collapse (mean) alphaZ e, by(j t teach_school)
	}
	if "`2'" == "preS" {
		keep if t==`3'
		collapse (mean) alphaZ e, by(j t teach_school_order)
	}


	tempfile tempjt
	save `tempjt', replace
	restore

	* choose sample

	if "`2'" == "loY" {
		keep if t!=`3'
	}
	if "`2'" == "preY" {
		keep if t<`3'
	}
	if "`2'" == "loS" {
		gen teach_school_t = teach_school if t==`3'
		bys j: egen teach_school_t_all = min(teach_school_t)
		keep if teach_school!=teach_school_t_all
	}
	if "`2'" == "preS" {
		gen teach_school_order_t = teach_school_order if t==`3'
		bys j: egen teach_school_order_t_all = min(teach_school_order_t)
		keep if teach_school_order<teach_school_order_t_all
	}
	if "`2'" == "wiS" {
		gen teach_school_t = teach_school if t==`3'
		bys j: egen teach_school_t_all = min(teach_school_t)
		keep if teach_school==teach_school_t_all
	}


	* construct some weights based on sample sizes

	gen n_ct1_n_ct2 = n_ct_m1*n_ct_m2
	gen n_ct1_sq = n_ct_m1^2
	gen n_ct2_sq = n_ct_m2^2

	gen Abar_m1_wgt = Abar_m1*n_ct_m1 // want to calculate mean residual, weighted by number of students
	gen Abar_m2_wgt = Abar_m2*n_ct_m2

	collapse (sum) n_ct_m1 n_ct_m2 n_ct1_n_ct2 n_ct1_sq n_ct2_sq Abar_m1_wgt Abar_m2_wgt (mean) alphaZ sigma* e, by(j t)

	gen Abar_m1 = Abar_m1_wgt / n_ct_m1
	gen Abar_m2 = Abar_m2_wgt / n_ct_m2

	gen w11 = n_ct1_sq / n_ct_m1^2
	gen w22 = n_ct2_sq / n_ct_m2^2
	gen w12 = n_ct1_n_ct2 / (n_ct_m1*n_ct_m2)

	assert w11<=1 if w11!=.
	assert w22<=1 if w22!=.
	assert w12<=1 if w12!=.

	* loop over teacher-years

	gen st = t-`3'
	gen abs_st = abs(st)
	gen drift_st = min(abs_st,`4')

	keep if abs_st<=`5'

	egen teacher_index = group(j)

	qui summ teacher_index
	local n = r(max)

	tempfile temploop tempstore

	save `temploop', replace

	local i = 1

	forv nn=1/`n' {

		use `temploop', clear

		keep if teacher_index==`nn'
		
		expand 2, gen(obsvar)
		gen m = obsvar+1
		
		gen Abar = (Abar_m1)*(m==1) + (Abar_m2)*(m==2)
		drop if Abar==.
		
		if _N > 0 {
		sort m t
		
		local nj = _N
		matrix A1 = J(`nj',1,.)
		matrix A2 = J(`nj',1,.)
		matrix gamma1 = J(`nj',1,.)
		matrix gamma2 = J(`nj',1,.)
		matrix Gamma1 = J(`nj',`nj',.)
		matrix Gamma2 = J(`nj',`nj',.)
		
		forv aa=1/`nj' {
			
			local m = m[`aa']
			local ta = t[`aa']
		
			local dd = drift_st[`aa']
			
			matrix A1[`aa',1] = Abar[`aa']
			matrix A2[`aa',1] = Abar[`aa']
		
			matrix gamma1[`aa',1] = sigma_m1m`m'_`dd'[`aa']
			matrix gamma2[`aa',1] = sigma_m`m'm2_`dd'[`aa']
		
			
			
			forv bb=1/`nj' {
				
				local tb = t[`bb']
				local mb = m[`bb']

				if `aa'==`bb' {
					matrix Gamma1[`aa',`aa'] = sigma2_muj`m'[`aa'] + sigma2_theta`m'[`aa'] * w`m'`m'[`aa'] + sigma2_eps`m'[`aa'] / n_ct_m`m'[`aa']
					matrix Gamma2[`aa',`aa'] = sigma2_muj`m'[`aa'] + sigma2_theta`m'[`aa'] * w`m'`m'[`aa'] + sigma2_eps`m'[`aa'] / n_ct_m`m'[`aa']
				}
				if `aa'!=`bb' & `ta'==`tb' {
					matrix Gamma1[`aa',`bb'] = sigma_muj1muj2[`aa'] + sigma_theta1theta2[`aa'] * w12[`aa']
					matrix Gamma2[`aa',`bb'] = sigma_muj1muj2[`aa'] + sigma_theta1theta2[`aa'] * w12[`aa']
				}
				
				if `aa'!=`bb' & `ta'!=`tb' {
		
					local drift_diff = min(abs(`ta'-`tb'),`4')
					matrix Gamma1[`aa',`bb'] = sigma_m1m`mb'_`drift_diff'[`aa']
					matrix Gamma2[`aa',`bb'] = sigma_m`mb'm2_`drift_diff'[`aa']
				
				}
			
			}
		}
		
		
		
		forv m=1/2 {
			gen mu_t_m`m'_hat = .
			if det(Gamma`m')>0 {
				matrix Gamma`m'_sym = 0.5 * (Gamma`m'+Gamma`m'')
				cap matrix psi = invsym(Gamma`m'_sym)*gamma`m'
				cap matrix VA = psi'*A`m'
				replace mu_t_m`m'_hat = VA[1,1]
			}
			
		}
		
		
		keep j mu_t_m1_hat mu_t_m2_hat
		keep if _n==1
		
		if `i'>1 {
			append using `tempstore'
		}
		sort j
		save `tempstore', replace
		local i = `i'+1
		}
		
	}

	use `tempstore', clear
	sort j

	merge 1:n j using `tempjt'
	drop if _m==2
	drop _m

	** add back in experience effect

	gen mu_jt_m1_hat = mu_t_m1_hat + alphaZ
	gen mu_jt_m2_hat = mu_t_m2_hat + alphaZ

	local extravars = "" 

	if "`2'" == "loS" | "`2'" == "wiS" {
		local extravars = "teach_school"
	}
	if "`2'" == "preS" {
		local extravars = "teach_school_order"
	}


	keep j t mu_jt_m1_hat mu_jt_m2_hat mu_t_m1_hat mu_t_m2_hat `extravars'
	sort j t

	ren (mu_jt_m1_hat mu_jt_m2_hat mu_t_m1_hat mu_t_m2_hat) ///
	 (mu_jt_m1_hat_`2' mu_jt_m2_hat_`2' mu_t_m1_hat_`2' mu_t_m2_hat_`2')


	save "$temp/blp_va_`1'_`2'_`3'_drift_behavioral", replace

	end

	* 1: subject
	* 2: sample description: "career" "loY" "loS" "preY" "preS" "wiS"
	* 3: sample specific value
	* 4: drift limit
	* 5: limit of sample used for measures

	foreach ss in ma { 
		
		use "$temp/VAstructural_`ss'_behavioral", clear

		qui summ t
		global tmin_`ss' = r(min)
		global tmax_`ss' = r(max)

		qui summ teach_school
		global smin_`ss' = r(min)
		global smax_`ss' = r(max)

		qui summ teach_school_order
		global sminorder_`ss' = r(min)
		global smaxorder_`ss' = r(max)
		
			local tmin = ${tmin_`ss'}
			local tmax = ${tmax_`ss'}
		
			forv tt=`tmin'/`tmax' {
	
				if `tt' > `tmin' {
					blp_va `ss' preY `tt' 8 6
				}

			}


	}

	
	** put data back together into one file

	foreach ss in ma {
		
			local tmin = ${tmin_`ss'}
			local tmax = ${tmax_`ss'}

			clear
			forv tt=`tmin'/`tmax' {
				if `tt' > `tmin' {
					append using "$temp/blp_va_`ss'_preY_`tt'_drift_B"
					replace t = `tt' if t==.
				}
			}
			save "$temp/blp_va_`ss'_preY_drift_B", replace

	}

	foreach ss in ma {
		use "$temp/VAschoolindices_`ss'_B", clear

		sort j t
		merge n:1 j t using "$temp/blp_va_`ss'_preY_drift_B"
		drop _m

		ren (*) =_`ss'
		ren (j_`ss' t_`ss' s_`ss') (j t s)
		
		save "$temp/blp_va_`ss'_drift_B", replace
	}

	foreach ss in ma {
		use "$temp/VAstructural_`ss'_B", clear
		gen Abar_m1_wgt = Abar_m1*n_ct_m1 // want to calculate mean residual, weighted by number of students
		gen Abar_m2_wgt = Abar_m2*n_ct_m2

		collapse (sum) n_ct_m1 n_ct_m2 Abar_m1_wgt Abar_m2_wgt (mean) sigma* e alphaZ vijt, by(j t s schlcode lea)

		gen Abar_m1 = Abar_m1_wgt / n_ct_m1
		gen Abar_m2 = Abar_m2_wgt / n_ct_m2
		
		drop *_wgt
		ren (Abar_m1 Abar_m2 n_ct_m1 n_ct_m2 sigma* e alphaZ vijt) =_`ss'
		
		tempfile temp`ss'
		save `temp`ss'', replace
	}


	use "$temp/blp_va_ma_drift_B"

	sort j t s
	merge 1:1 j t s using `tempma'
	drop _m

	foreach ss in ma {
		foreach nn in preY {
			gen mu_jt_hat_`nn'_`ss' = mu_jt_m1_hat_`nn'_`ss' * p_m1_`ss' +  mu_jt_m2_hat_`nn'_`ss' * p_m2_`ss'
		}
		gen Abar_`ss' = Abar_m1_`ss' * p_m1_`ss' + Abar_m2_`ss' * p_m2_`ss'
	}

	save "$basedata/va_estimates_drift_behavioral", replace
}
